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Abstract 

q ■ We present a new module of the micrOMEGAs package for the calculation of 

QJ . WIMP-nuclei elastic scattering cross sections relevant for the direct detection of 

Q \ dark matter through its interaction with nuclei in a large detector. With this new 

' module, the computation of the direct detection rate is performed automatically 

for a generic model of new physics which contains a WIMP candidate. This model 
needs to be implemented within micrOMEGAs 2 . 2. 

Ph; 1 Introduction 

The existence of an important cold dark matter (CDM) component has been firmly es- 
! tablished by cosmological observations in the last few years notably by SDSS P] and 
WMAP [2]. A leading candidate for CDM is a new weakly interacting massive particle 
(WIMP). This WIMP must be stable. Such particles arise naturally in many extensions 
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of the standard model [3] from the minimal supersymmetric standard model [H [S] to 
models of extra dimensions P, El [9], little Higgs models [10] or models with extended 
gauge [nj 02] or Higgs sectors [131 E] • 111 these models the dark matter (DM) candidate 
can be either a Majorana fermion, a Dirac fermion, a vector boson or a scalar. Their 
masses range anywhere from a few GeV's to a few TeV's. 

Astroparticle experiments are actively pursuing searches for WIMP DM candidates 
either directly through detection of elastic scattering of the WIMP with the nuclei in a 
large detector or indirectly trough detection of products of DM annihilation (photons, 
positrons, neutrinos or antiprotons) in the Galaxy or in the Sun. 

In direct detection, one measures the recoil energy deposited by the scattering of 
WIMPs(x)0 with the nuclei. Generically WIMP-nuclei interactions can be split into spin 
independent (scalar) and spin dependent interations. The scalar interactions add co- 
herently in the nucleus so heavy nuclei offer the best sensitivity. On the other hand, 
spin dependent interactions rely mainly on one unpaired nucleon and therefore dominate 
over scalar interactions only for light nuclei unless scalar interactions are themselves sup- 
pressed. In both cases, the cross-section for the WIMP nuclei interaction is typically 
low, so large detectors are required. Many experiments involving a variety of nuclei have 
been set up or are being planned. Detectors made of heavy nuclei (for example Germa- 
nium or Xenon) currently in operation include Edelweiss [15], DAMA |16| . CDMS [T7] . 



1 Hcre we use x to designate the DM candidate whether a fermion. scalar or vector boson. 
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Xenon [IB], Zeplin [20], Warp [21] and KIMS [22]- Upgrades and new projects such as 
Genius, Xmass |23j, CLEAN [21], ArDM [25] and Eureca [26] have been proposed as 
well. Detectors made of light nuclei which are sensitive mainly to the spin dependent 
interaction include Simple [22], Picasso [28], Tokyo/NaF [29] and NAIAD [3D]. The 
latter having one light (Na) and one heavy (I) target nuclei is actually sensitive to both 
spin-dependent (SD) and spin-independent (SI) interactions. Larger versions of existing de- 
tectors and new projects are also proposed, for example those operating with 3 He [32l [33] . 
Note that heavy nuclei although best for probing the scalar interaction have also a sen- 
sitivity to spin dependent interactions because of their odd-A isotopes. Currently the 
sensitivity of both types of detectors for spin dependent interactions is similar. Further- 
more, different nuclei offer a sensitivity to spin dependent interactions on protons (for 
odd-proton nuclei such as 23 Na, 127 I or 19 F) or neutrons (for odd-neutron nuclei such as 
or 29 Si 73 Ge, 129 Xe). 

Only one experiment, DAMA, has reported a positive signal consistent with an anni- 
hilation cross-section a xn « 0.2 — 1. x 10~ 5 pb for a WIMP mass around 30-100 GeV [34] . 
Other experiments, such as Edelweiss, CDMS, or Xenon have only set an upper limit on 
the WIMP-nucleon annihilation cross-section_|. The best limits were reported recently 
by Xenon, erg « 4 x 10^ 8 pb for a WIMP mass around 30 GeV [32] and by CDMS, 
erg ~ 4.6 x 10~ 8 pb for a WIMP mass around 60 GeV [19]. These values already probe 
a fraction of the parameter space of the most popular CDM candidate, the constrained 
minimal supersymmetric standard model (CMSSM) [3B1 EZ] or some of its extensions [TT] 
and poses severe constraints on a model with a Dirac right-handed neutrino [SSj- The 
search for WIMPs will continue with larger detectors (around 100kg) planning to reach a 
level of erg ~ 10~ 9 pb by 2010. By 2015, improved large detectors, around 1 ton, should 
go below the 10~ 10 pb level, for example Warp, Xenon, Eureca [26] or SuperCDMS [39] . 
For spin dependent interactions, the limits for neutrons from Zeplin [H] and CDMS [12] 
have recently been superseded by Xenon [4"U] , a^ D ~ 5. x 10~ 3 pb, while for protons 
the best limit from direct detection experiments was set by KIMS, of D « 0.18pb [31]. 
Indirect detection experiments looking for an excess of muon neutrinos from WIMP anni- 
hilations such as Super-Kamiokande have also set a stringent limit on the WIMP proton 
cross section, <rf D « 3. X 10~ 3 pb [12] • These limits are at the level or below the positive 
signal reported by DAMA |34j . 

The calculation of the cross-section for WIMP scattering on a nucleon have been ob- 
tained at tree-level for different DM candidates: neutralinos in supersymmetry(for reviews 
see [B],[l5]), gauge bosons in UED models [8] or in little Higgs models [IB], right-handed 
neutrinos [17] or scalars [131 113 US]. Implications of the direct detection experiments 
on DM models have been explored for quite some time [SUl EH E21 E31 EH ESI EB]- In 
the MSSM, the most complete calculation of the neutralino nucleon scattering is the one 
of Drees and Nojiri [57] that includes higher-order effects from twist- 2 operators. Public 
codes for DM in the MSSM such as DarkSUSY [5H] and Isajet [50] both follow this ap- 
proach for calculating neutralino nucleon scattering. On the other hand, micrOMEGAs 2 . 2, 
a code primarily designed for the calculation of DM relic abundance did not, up to now, 
provide a module for the computation of the direct detection rate even though an im- 
proved tree-level computation within the MSSM has been performed for some time [59] . 
This is the gap we intend to fill. Here we describe the implementation of the direct DM 

2 Possibilities for reconciling DAMA results with other experiments have been considered in Ref. [35] 
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detection rate within micrOMEGAs 2.2. 

Many ingredients enter the calculation of the direct detection rate and cover both 
astroparticle, particle and nuclear physics aspects. We need to know the WIMP density 
and the velocity distribution near the Earth. Since the WIMP have small velocities, it 
means that the momentum transfer, Q 2 , is very small as compared to the masses of the 
WIMP and/or nuclei. The detection rate depends of course on the WIMP nucleus cross 
section. To arrive at the x-nucleus cross section one has to first compute the interaction 
at the more fundamental level, that is at the quark level. The different matrix elements 
for \q interactions that capture the dynamics of the model in a perturbative way have 
to be converted into effective couplings of WIMPs to protons and nucleons. Finally we 
have to sum the proton and neutron contribution and turn this into a cross section at 
the nuclear level. The recoil spectrum of the nuclei depends on the velocity distribution 
and, in view of the low Q 2 , is contained in the elastic form factor of the nucleus. In this 
manual, we describe all these steps. Even though many of these steps are not new, it is 
necessary to understand how all the pieces are implemented in the code. 

An important point to emphasize is that micrOMEGAs, contrary to other public codes, 
is not restricted to the supersymmetric model with a neutralino DM but is applicable 
to a generic model of new physics for DM E). We have already in previous versions set 
up the code so that any model can be implemented to give the relic density of DM, the 
indirect detection rate and cross-sections relevant for collider applications. In the same 
spirit, the calculation of direct DM detection rate in nuclei is also performed for generic 
models of DM. More precisely the tree-level cross-section is computed in any model and 
dominant QCD corrections are taken into account. Other higher order corrections such as 
the threshold corrections to Higgs quark vertices are model dependent and are provided 
only for the MSSM and its extensions (CPVMSSM, NMSSM). The steps that go from the 
automatic computation of the cross-section for WIMP scattering on quarks to a detection 
rate in a large detector follow standard approaches [Hj. In the spirit of the modular 
approach of micrOMEGAs, different nuclear form factors or WIMPs velocity distribution 
can easily be implemented by the user. 

The paper is organised as follows, we first review the computation of the scattering 
rate for DM on a point-like nucleus starting from an effective Lagrangian for nucleon- 
WIMP interactions. We then show in Section 3 how to relate these to the quark- WIMP 
interactions and describe the method used to reconstruct the effective Lagrangian for both 
scalar and spin dependent interactions. We also describe the treatment of dominant QCD 
corrections. The computation of the recoil distribution for WIMPs scattering on nuclei 
taking into account nuclear form factors and velocity distribution of WIMPs follows in 
Section 4. The functions available in this new module are described in Section 5. Section 
6 is devoted to sample results and comparisons with other codes. The treatment of box 
diagrams is described in Appendix A and some details about nucleus spin dependent form 
factors are gathered in Appendix B. 

3 A DM code, DM++, that can also be applied to any new physics model has been developed by |3B] . 
This code is not yet publicly available. 
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2 Elastic scattering of WIMPs on point-like nuclei 



The standard formalism for evaluating WIMP nuclei cross sections was reviewed in [Sj 
with special emphasis on the case of the neutralino DM. Here we first describe the cal- 
culation of DM elastic scattering on point-like nuclei taking the Majorana fermion as an 
example then consider different types of DM candidates. The relation with the scattering 
rates on quarks, computed automatically in our code will be presented in the next section. 

The velocity of DM particles near the Earth should be of the same order as the orbital 
velocity of the Sun, v ~ 0.001c. Because of this small velocity, the momentum transfer 
is very small as compared to the masses of the WIMP and/or nuclei. For example for 
typical masses of WIMP, M x m 100 GeV and of nuclei, Ma ~ 100 GeV, the maximum 
transfer momentum is 

JTZqi = 2v ' ] , lx ' ]l :\ « lOOMeV « 0.5/m- 1 . (1) 
v ^ M X + M A J w 

Thus all WIMP nucleon elastic cross sections can be calculated in the limit of zero mo- 
mentum transfer. The cross sections for scattering on nuclei are obtained from the WIMP 
nucleon cross sections after folding in the nuclei form factors, these form factors depend 
on the momentum transfer. 

2.1 Scattering rate on point-like nuclei- the case of a Majorana 
fermion 

In the non-relativistic limit, WIMP-nucleon elastic amplitudes can be divided into two 
classes, the scalar or spin independent interaction and the axial- vector or spin dependent 
interaction. For a spin 1/2 nucleon, interactions corresponding to multipole will clearly 
vanish in the zero momentum limit. In the familiar case of a Majorana fermion, the 
effective Lagrangian reads [60J 



+ ^xl^lb^x^Nl^N + iN^xliil^x^Nl^l^N (2) 

In the zero momentum transfer limit the operator ipj^ip vanishes while only the space 
component of V'TsT/^ an d the time component of ip'j^ remain. Thus the operators 
are suppressed in the limit of small momentum transfer by factors of order q 2 /m 2 N and/or 
q 2 /m x where q 2 = —Q 2 . We will ignore these operators. Therefore only one operator 
survives in each class, Xn for SI and £/v for SD. Of course in a specific model it is possible 
that the coefficient Xn is much smaller than one of the coefficients Kj, in which case the 
operator may contribute at the same level as X^ despite the Q 2 suppression. However 
in this case we expect very small rates, much below the experimental sensitivities. 

Spin independent interactions 

For SI interactions with nucleons the effective Lagrangian thus reads 

C SI = X N ipijj x ip N ip N (3) 
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where N = p,n. The squared amplitude for a nucleon after averaging (summing) over 
the polarization of incoming (outgoing) particles is, 

|^ / | 2 = 64(A, V M X M 7V ) 2 (4) 

where is the nucleon mass. Scalar and vector WIMP-nucleon interactions naturally 
induce scalar and vector WIMP-nuclei interactions. Summing on proton and neutron 
amplitudes gives for WIMP-nucleus interaction at rest, 

|^| 2 = 64M 2 M 2 (A P Z + X n (A - Z)) 2 (5) 

where Z is the nucleus charge and A the total number of nucleons. It leads to the cross 
section for a WIMP scattering at rest from a point-like nucleus 

4u 2 

oSi = J*(\ pZ + X n (A-Z)) 2 (6) 

7T 

where \i x is the WIMP-nucleus reduced mass, fi x = M x Ma/(M x + M A ). Note that the 
nucleon cross-section adds coherently so that there is a strong enhancement for large 
nuclei, oc A 2 when \, p « A n . 

Spin dependent interactions 

The effective Lagrangian for spin dependent interactions of a Majorana fermion at 
zero momentum transfer reads 

£ SD = Zn^ x 1s1^ n 1^n (7) 

It leads to the squared amplitude 

\A 6 N D \ 2 = 192(£ N M X M N ) 2 (8) 

In order to get the amplitudes for nuclei we have to sum the spin currents produced 
by the protons and neutrons separately. Since we know that the spins of protons with 
the same orbital state should be opposite, we expect strong compensation of currents 
produced by protons as well as those produced by neutrons. First note that for interactions 
at rest, the 70 component of the pseudovector current, Eq. [7J vanishes. The resulting 
interaction i^^^iip leads to a three dimensional vector current. This vector current has 
to be proportional to the angular momentum J. We can write for nuclei 

J$ = S%f A /\J A \ (9) 

where Sfi are the expectation value of the spin content of the nucleon N in a nucleus with 
A nucleons. By definition, for protons and neutrons S£ = S™ = 0.5 and S™ = = 0. 

The second peculiarity of the SD case is a non-trivial summation over spins. Because 
the matrix element is proportional to J A the summation over spin states in a nucleus gives 
a factor 

Yl Yl Yl < s x\4K >< s x\ J x\ s x >< s a\J a W a >< s' A \J l A \s A > 
= E M^4M j W) = (2J x + 1)^x + 1)(2^i+1)Ja(Ja + 1)/3 (10) 

l<fc,K3 
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here we use s, s' for labelling polarization states and J a refers to the angular momentum 
of a nucleus with A nucleons. After averaging over initial polarizations, a factor (2J X + 
1)(2Ja + 1) will cancel out. 

Taking into account the spin currents structure (jHj) and the J dependence ( TTOT) we can 
write the WIMP-nucleus squared amplitudes as 

\A SD \ 2 = 256^±I foSf + ^) 2 M 2 M 2 A (11) 

J A 

This reduces to Eq. [8] in the special case of the nucleon and leads to the cross section at 
rest for a point-like nucleus [H], 

\tpSf + tnSff (12) 



16/4 J A + 1 



SD 
'0 



7T X) 



The quantities are obtained from nuclear calculations or from simple nuclear mod- 
els, such as the odd-group model. They are estimated to be « 0.5 for a nuclei with an odd 
number of protons or neutrons and ~ for an even number. Thus no strong enhancement 
is expected for SD interactions in nuclei. The treatment of the nuclei form factors taking 
into account the momentum dependence will be discussed in section 5.2. 



2.2 Generalization to other DM candidates 

To derive the formulae for elastic scattering on point-like nuclei, we started from the ef- 
fective WIMP-nucleon Lagrangian ( |3|71) written for a Majorana WIMP. In fact these can 
be generalized to all types of WIMPs . Our aim it to give the generic form of the effec- 
tive Lagrangians for a fermionic, scalar and vectorial WIMP including the possibility of 
complex fields. In all cases we define the effective Lagrangian such that the normalization 
conditions ( 141IB1) are satisfied. Here we write only operators that contribute at q 2 = 0. 
As we have argued for Majorana fermions, other operators are suppressed by q 2 / m x (A) 
and can potentially be of the same order as the operators we consider only when both 
contributions to the scattering cross section are small. Note that in the case of a complex 
field, x and X have in general different cross-sections. For each type of interaction, SI (SD) 
one can then construct two operators, one that is even with respect to x ~ X interchange, 
^N,e(£,N,e) and is the only remaining operator for Majorana's and another that is odd, 
^N,o (£n,o)- 

For a fermion field the most general Lagrangian reads 

+ 6v,e'0x757M^X^iV757 / ViV ~ ^N^x^^X^N^IpN (13) 

This Lagrangian leads to matrix elements which satisfy the normalization conditions 
Eq. KED with 

Atv = and £jv = ^ (14) 

where the +(— ) signs correspond to WIMP (anti-WIMP) interactions. The special case of 
a self-conjugated WIMP such as the Majorana fermion is recovered when A^v o = 6v o — ► 
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and A at = Ajv e , £n = £jv e - Note the factor 2 difference between the operator for Majorana 
and Dirac fermion field, compare Eq. [3] and Eq. [T31 Note that the antisymmetric tensor 
a^v current effectively reduces to a vector interaction since in the non-relativistic limit 
only the spatial components contribute. 

For a scalar field only SI interactions are possible, for the general case of a complex 
scalar, 

C s = 2\ NtE M x (p x (j)* x ij N ^ N + i\N,o{d^<p x 4> x - (^xd^D^Nl^N (15) 

The squared amplitude is normalized as for Majorana fermions, (j3J), with the condition 
(j!4p for complex scalars. Again the case of the real scalar corresponds to Atv, = and 
A at = Ajv,e- Note that the four- dimensional vector current ip^j^ipN actually leads to a 
scalar interaction because only the zeroth component of this current does not vanish in 
the non-relativistic limit. 

Finally for a complex vector field, 

C v = 2X Nie M x A w A x ^p N ip N + X N A A T 9 ^ A x,c - A x a d^A xa ]i} Nl ^ N 
+ V6Ue(d a A xf3 A x ^ - A* xp d a A xl )e a ^ Nl ^ N 

+ i ^ Nt0 (A x ,A* xu - A Xfl A xv )^ N a^ N (16) 

with in the special case of a real vector field Atv, = 6v,o = 0. Again the couplings are 
normalized as for the fermion case both for real (JU [H]) and complex fields (j!4p . 

In micrOMEGAs 2 . 2 we assume that DM particles and anti-particles have the same 
density. We do not consider the case where this symmetry is broken by CP violation. 
Under this assumption, the event rate for WIMP scattering in a large detector is obtained 
after averaging over \- and x-nucleus cross-sections. 

3 WIMP elastic scattering on quarks 

The matrix elements for WIMP nucleon interactions are related to the more fundamental 
matrix elements for WIMP quarks interactions. These matrix elements can be easily 
written in a given model. To handle a generic model we rather expand WIMP quark 
interactions over a set of basic point-like operators. Only a few operators are necessary 
in the q 2 — > limit. 

The operators that are non-zero in the non-relativistic limit are similar to the operators 
introduced for nucleons in Eq. [T3lll5|16l with ip^ —>■ ifj q . Those operators are listed in 
Table[T]for either spin-independent or spin-dependent interactions of a scalar (0 X ) fermion 
(\l/ x ) or vector (Alt) WIMP. Note that for the latter we use the unitary gauge. Of course 
a scalar WIMP can only contribute to spin independent interactions. As for nucleons, 
the operators are separated in two classes - odd and even, depending on the symmetry 
properties with respect to quark - anti-quark exchange. For real WIMPs, odd operators 
are zero by definition. 

The operators in Tables [1] have the same normalization as the operators for WIMP 
nucleon interactions, that is the squared matrix element for q\ — > qx SI interactions after 
summing/averaging over polarizations of outgoing/incoming particles are 

\A?\* = <* ^ + 2 ^ o) M x m q y (17) 
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Table 1: Operators for WIMP - quark interactions. 





WIMP 


Even 


Odd 




Spin 


operators 


operators 


SI 




1/2 
1 


2MJ> x $$ q ^ q 


+i\ q , (AT 9 » A x,« ~ A x a d,A xa )^ ql ^ g 


SD 


1/2 
1 


V6(d a A* x pA xu — A* x pd a A xv ) 


i~(A xtl A* xu — A xll A xu )ip 'ip q 



where A g e (Ag ) are the coefficients of the even (odd) operators in Table 1. In the special 
case of a pure neutral WIMP, {\, e + A 9)0 )/2 — > X q . Similarly for SD interactions the 
normalization is specified in Eq. [HJ 

3.1 Numerical approach to operator expansion 

The micrOMEGAs 2 . 2 package contains an automatic generator of matrix elements, CalcHEP [76] . 
For each model, micrOMEGAs has a complete list of Feynman rules which specify the model 
in the CalcHEP format. Therefore CalcHEP can generate automatically amplitudes for 
XQ ~ > XQ f° r an Y kinematics. Here we need to generate these matrix elements at low 
energy in a format that can be easily and automatically turned into effective operators 
for x-nucleon interactions. So we want to expand automatically the Lagrangian of the 
model in terms of local operators and extract the coefficients of the low energy effective 
WIMP quark Lagrangian 

£eff(x) = Ks6 q>s {x) + Us6' q ,s{x) (18) 
q,s 

where q is a label for quark, s is a label for even and odd operators and O q)S (0' q , s ) 
are the spin independent (dependent) operators in Table [U Traditionally the coefficients 
\,sj £q,s are evaluated symbolically using Fiertz identities in the limit q 2 <S M\. Instead in 
micrOMEGAs 2.2 these coefficients are estimated numerically using projection operators. 
First we compute special matrix elements for \q — > \q scattering at zero momentum 
transfer. For this we need to add to the model Lagrangian the projection operators defined 
in Table 1. The interference between one projection operator and the effective vertices will 
single out either the spin dependent or spin independent contribution, since the effective 
Lagrangian is written in an orthogonal basis. To further separate the coefficient of the 
even and odd operators, we compute both \Q ~^ XQ an d XQ ~^ XQ matrix elements. We 
use the fact that for a given projection operator the interference term in the squared 
matrix elements for quarks are identical to those for antiquarks for even operators and 



S 



have opposite signs for odd operators. Thus, taking into account the normalisation, 

(?(pi), x(P2)|O g , e d gi e|g(pi), x(P2)) 
A _ A _ -»(g(pi), x(P2)|^C g , e |g(pi), x(P2)) / 19 n 

(?(Pl), x(P2)|O g ,eO,, B |g(pi), x(P2)) 

where the S'-matrix, S = 1 — iC is obtained from the complete Lagrangian at the quark 
level. 

To implement this procedure, a new model file is created automatically in micrOMEGAs 2 . 
This file contains the model file of a given model as well as the auxiliary vertices of Ta- 
ble [U There is no need for the user to implement these auxiliary vertices. CalcHEP then 
generates and calculates symbolically all diagrams for WIMP - quark/anti-quark elastic 
scattering keeping only the squared diagrams which contains one normal vertex and one 
auxiliary vertex. This corresponds to the matrix elements in Eq. [TH1 Note that in the file 
that defines the model all quarks should be defined as massive particles. Vertices that 
depend on light quark masses cannot be neglected, for example the couplings of Higgs to 
light quarks, since the dominant term for WIMP quark scalar interactions is proportional 
to quark masses. In particular in SUSY models masses of the first and second generation 
fermions must be included. When converting to WIMP nucleon interactions this quark 
mass will get replaced by a nucleon mass so that the nucleon cross section is independent 
of the light quark masses, see Section 13.31 Because the amplitude for WIMP scattering 
on light quarks is proportionnal to a small quark masses, one has to be wary of numerical 
instabilities. To avoid these we consider the amplitude for q(pi)x(P2) ~> q(P3)x(Pi) an d 
write the squared matrix element is terms of the dot products p\.p2 and p\.p%. The ampli- 
tude depends explicitly on the small momenta p\ thus avoiding the numerical instabilities 
in matrix elements that were found with other kinematics. 



3.2 Contribution of tensor (twist-2) operators. 

There are other point-like operators which produce the same SI and SD amplitudes at zero 
momentum. They are the operators containing field derivatives. Here we consider the 
contributions of the symmetric traceless tensor operator to the WIMP nucleon scattering. 
A complete treatment of twist-2 operators in neutralino nucleon elastic scattering in 
the MSSM was first presented in [57]. In the MSSM tensor operators come from the 
momentum expansion of the denominator in the squark exchange diagram and are usually 
suppressed by a factor M N /(Mg — M x ) with respect to the main contribution. Thus they 
are expected to be mainly relevant when the squark is not much heavier than the WIMP, 
for example in the coannihilation region. 

In the MSSM, the tensor operator reads [57] . 

O q ,t = ^(X7/AX)<T (20) 

£>f = q(^~d v - -f*d v + 7""0" - -f^d* + im q g^)q (21) 

where 0% u is a twist-2 operator (recall that twist is defined as dimension — spin). Note 
that this tensor operator is generically found in other models where new coloured particles 
couple to WIMP and quarks. 
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Figure 1: Diagrams that contribute to WIMP-gluon interaction via quark loops in the 
MSSM. 

To extract automatically the coefficient of the tensor operator in the low energy WIMP 
quark Lagrangian we consider forward scattering at small momentum. Indeed for such 
kinematics the matrix element of the product of two scalar operators does not depend on 
the collision while the product of scalar and tensor does since 

< g(Pi),x(P2)|0,,A,e|g(pi),x(P2) >= -32m q M x (4( Pl .p 2 ) 2 - m 2 q M 2 x ) (22) 

Thus the tensor operator can be extracted by evaluating numerically the second derivative 
of 

< q{pi),x{P2)\SO q MPi),x{P2) > (23) 

with respect to p\.pi- Note that this trick does not require that one implements the 
tensor operators in the CalcHEP model. Since the tensor operators contribute also to 
the amplitude at rest, after extracting the coefficient of the tensor operator it should be 
subtracted from Eq. [19] to isolate the coefficient of the scalar operator. 

In typical MSSM models the correction due to an accurate treatment of twist-2 op- 
erators is less than 1%, in special cases though, for example for a small mass difference 
between the squark and the WIMP the contribution can be larger. The contribution 
of higher spin twist-2 operators are in general further suppressed and are not included 
here. In the MSSM we have compared numerically the cross sections obtained with our 
method for extracting the effective operators, including the twist-2 operators, with the 
complete results of Drees and Nojiri [57]. For this we have implemented their analytical 
results, correcting a sign in the SD amplitude between Z and squark exchange. We found 
numerical agreement at the percent level. 

3.3 Nucleon form factors 

In order to convert WIMP-quark amplitudes to WIMP-nucleon amplitudes, we need to 
know the values of the quark currents inside the nucleon. The operators for these quark 
currents in a nucleon can be extracted from experiment or estimated theoretically. In the 
following we give estimates for the coefficients of the four quark currents listed in Table Q] 

3.3.1 Scalar coefficients 

The scalar ip q ip q current characteristic of the Si-even interaction depends on the total num- 
ber of quarks and anti-quarks in the nucleon. The operator {N\m q ip ^ q \N) is interpreted 



10 



as the contribution of quark q to the nucleon mass, M^, 



(N\m q ^ q \N) = f^M N (24) 
where the coefficients f q relate nucleon and quark operators, 

9=1,6 

Note there is no explicit dependence on the quark mass in the cross section for WIMP 
nucleon scattering. Indeed the quark mass term gets transformed into a nucleon mass. 
For heavy quarks, Q, the parameter fg is induced via gluon exchange with the nucleon, 
see section 13751 and [77] 

(26) 

The coefficients f q can be determined using the value of the light quark masses extracted 
from baryon masses, the ratio of the quantities B q = (N\qq\N) for u, d and s quarks 
and from the value of the pion-nucleon sigma-term. It is the latter that has the largest 
uncertainty, see for example [TJJ EQ]. For light quark masses ratio we take 

Tfl TfX 

— = 0.553 ± 0.043, — = 18.9 ±0.8 (27) 
m d m d 




We define the quantities 



|Hr«l-49 and y = (28) 



where y denotes the fractional strange quark content of the nucleon. y is obtained from 
the pion-nucleon sigma term and the quantity ctq which is related to the size of the SU(3) 
symmetry breaking effect [81] 



tr^ = - \ d ' (B u + B d ); a = \ d - (B u + B d - 2B S ) (29) 



which implies 



1 - (30) 

O'nN 



Recent analyses suggest that 

a vN = 55 - 73 MeV and a = 35 ± 5 MeV (31) 



where o"o is estimated from chiral perturbation theory [81] or from baryon mass differ- 
ences [83]. Defining 

B, 2z -{z- l)y 
a =g- d = 2 + (.-!)» (32) 



the parameters for light quarks in the proton write 



Jd- (1 + mh+aV Tu ~ m„ aJd ' Is ~ (I + ^)m n m d [ ' 
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and in the neutron 

2a nN a fn ^rnul fn fn = <t*nV rn s 
Jd (l + ^) mn (l + a )' J - md a Jd > Js (l + ^)m n m d 1 } 

As default values we take <jq = 35 MeV and ovjv = 55 MeV which lead to 

f d = 0.033, fl = 0.023, /f = 0.26 

f2 = 0.042, / u n = 0.018, /r = 0.26 (35) 

In micrOMEGAs 2.2, the values for these coefficients can be changed directly or through 
modification of oo and a n N and the quark mass ratios. In all cases the parameters for 
heavy quarks will be recalculated using Eq. [261 Note that is typically larger than 
the value used in earlier analyses, /| = 0.118 — 0.14. This is mainly due to an increase 
in a n N which was centered around 45 MeV [HH ES]- This large correction to can 
lead to an increase by a factor 2-6 in the spin independent cross-section for nucleons [51] . 
Furthermore even with the new estimate of a n N, large uncertainties remain. 



3.3.2 Vector coefficients 

The vector if) jy^ipq current in the Si-odd interaction is responsible for the difference be- 
tween xiV and xN cross sections. The interpretation of this current is very simple. It 
counts the number of quarks minus the number of anti-quarks in the nucleon, that is the 
number of valence quarks. This current is the only one that does not suffer from theoret- 
ical uncertainties when going from the WIMP- quark interaction to the WIMP-nucleon 
interaction. Indeed only valence quarks contribute to the vector current so that 

q=u,d 

with f Vu =2,f Vd = !,/»=!,/« =2. 



3.3.3 Axial-vector coefficients 

The axial- vector current ^ q ^^^ q is responsible for spin dependent interactions. It counts 
the total spin of quarks and anti-quarks q in the nucleon. Operators for axial-vector 
interactions in the nucleon are related to those involving quarks, 

6v, e = (37) 

q=u,d,s 

with 

2 8(t Aq N = (N® q w<J> q \N) (38) 

Here s M is the nucleon spin and Aq N are extracted from lepton-proton scattering data. 
The strange contribution to the spin of the nucleon, as measured by EMC and SMC 
turned out to be much larger than expected from the naive quark model [EE]. This leads 
to the following estimates for the light quark contributions in the proton which have been 
used in many analyses of DM spin dependent interactions, 

A£ = 0.78 ±0.02, = -0.48 ±0.02, = -0.15 ± 0.02 (39) 



12 



These early results have qualitatively been confirmed by HERMES [HZ] and also by COM- 
PASS [88] . As default values for the axial- vector coefficients we use the latest determina- 
tion of the light quark contributions [87J, 

A p = 0.842 ±0.012, A^ = -0.427 ±0.013, A^ = -0.085 ± 0.018 (40) 

These results are obtained in the limit of SU(3)f symmetry. It is argued that flavour sym- 
metry breaking effects might lead to an additional shift in the strange quark contribution, 
thus making this value compatible with the value extracted from EMC and SMC [86] . 
We neglect the heavy quark contribution to the proton spin since they have been shown 
to be small [89] . The neutron quantities are simply obtained by an isospin rotation 

A n u = A p , A2 = A£, A: = A p (41) 

Note that because there can be a cancellation between the Ag's when summing over 
light quarks in the nucleon, there can be a strong reduction in the coupling of WIMPs to 
neutrons or protons when varying the quark coefficients within the error bars. 

3.3.4 Coefficients of the a^ u term 

The tensor current, tj; a^ipq in SD-odd interactions is responsible for the difference be- 
tween \N and x.N spin-dependent cross sections. This current can be interpreted as the 
difference between the spin of quarks and the spin of anti-quarks in nucleons. Recent 
measurements by COMPASS [9DJ and HERMES [9T] indicate that the antiquark con- 
tribution (Su N + 5d N ) is compatible with zero. Forthcoming COMPASS data on muon 
scattering off a proton target will provide a separate determination of Au p and Ad p [90] . 
If confirmed to be zero, it would mean that the coefficients for the tensor interaction are 
identical to those for the axial vector interaction. The tensor coefficients have also been 
computed on the lattice, for example Ref. [92] gives 

5 P U = 0.84, 5 p d = -0.23, 5 P = -0.05. (42) 

These results were later confirmed by another lattice calculation of the LHPC collabo- 
ration [93]. We do not attempt to estimate the uncertainty associated with the lattice 
computation of the tensor coefficients. Note that both collaborations also compute the 
axial-vector coefficients on the lattice and they find values compatible (with larger un- 
certainties) to the one obtained from COMPASS and HERMES data. In our code it 
is possible to modify the input values for the tensor coefficients independently of the 
axial-vector coefficients. 

3.4 Twist-2 coefficients 

The nucleon form factors for twist-2 operators are well known in QCD. For the operator 
0%, eq. [2H we have [HZ] 

< N(p)\0%\N(p) >= (p»p v /M N - g^M N /A) [ (q(x) + q{x))xdx (43) 

Jo 

where q(x) and q(x) are parton distribution functions calculated at the scale Q = Mg—M x . 
We use the CTEQ6L distribution functions [93]. The computation of the box diagrams 
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for neutralino gluon scattering in [57] leads to a loop improved twist-2 tensor form 
factor containing gluon densities instead of quark parton densities. However these authors 
argued that the tree-level approach was more robust for c and b quarks because of large 
logs in the loop result, log(Q/m q ). In our package we use the tree level formulae for both 
c and b quarks and we neglect the twist-2 form factor for t quarks. @ 

3.5 Gluon contribution and QCD corrections to the scalar am- 
plitude 

The nucleon consists of light quarks and gluons, nevertheless one can also consider in- 
teractions of WIMPs with heavy quarks inside the nucleon. These effectively come into 
play since heavy quark loops contribute to the interactions of WIMP with gluons, see for 
example the diagrams in Fig. [TJ It is well known that when Higgs exchange dominates, 
see for example [57], WIMP gluon interactions via heavy quark loops can be taken into 
account by considering WIMP interactions with heavy quarks together with an estimation 
of the heavy quark condensates in nucleons. Furthermore dominant QCD corrections can 
also be taken into account in that case [95] . 

The anomaly of the trace of energy-momentum tensor in QCD implies [77] 

M N (N\N) = (N\ m^M 1 + 7) + (t^HG^GHAO (44) 

q<nt s 

where 7 is the anomalous dimension of the quark field operator, a s the strong coupling 
constant, the gluon field tensor and (3 n f = — a^/47r(ll — 2rij/3+a s /47r(102— 38nj/3)). 
In the leading order approximation for three flavours, Eq. (1441) is simplified to 

M N (N\N) = (N\ m^^-^afi^G^N) (45) 

q=u,d,s 

Comparing Eq. HH for n f and rif + 1 one finds the contribution of one heavy quark 
flavour to the nucleon mass, relating the heavy quark content of the nucleon to the gluon 
condensate 

(N\m Q ^ Q \N) = - *f (N\a.G lw Gr\N) (46) 

2a;(l + 7) 

where A/3 is a contribution of one quark flavor to the QCD /3-function. This formula 
agrees with the effective HG^ V G^ V vertex at small q 2 [96J. Up to order a 2 ,, 

(N\m Q ^ Q \N) = + Ua j mQ) ){N\a s G, v GnN) (47) 

Keeping only the leading order and combining with Eq. |45j will lead to the usual relation 
Eq. [2S Note that the NLO terms in Eq. [H partially cancel the effect of the NLO 
corrections in Eq. H7J so that the QCD corrections to Eq. [2S are small. See also [HTJ for 
an alternative estimate of the heavy quark content of the nucleon. 

4 Note that we evaluate the parton distributions at a scale Q = Mq—M x rather than Q = (M~— M^) 1 / 2 
as suggested in [57]. In this case we find that the loop and tree-level approaches agree reasonably well 
for 6-squarks when — M x m 20 GeV. Larger mass differences give a very small contribution from the 
twist-2 term while smaller mass differences typically lead to a low value for the DM relic density. 
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The simplest way to take into account dominant QCD corrections to Higgs exchange is 
then to consider WIMP heavy quark interactions through Higgs exchange and introduce 
an effective vertex for heavy quarks in the nucleon with Eq. [26] modified to include one- 
loop QCD corrections, Eq. H7J The equivalence of this approach with the description of 
the Higgs coupling to the nucleon through gluons is confirmed by a direct computation 
of the triangle diagram of Figfj] in the limit where Q 2 <C Mq. Recall that the typical 
transfer momentum is Q « 100 MeV, Eq. [TJ Note that for light quarks the corrections 
that would arise from their contribution to the triangle diagram that couples a Higgs to 
gluon are all absorbed into the definition of the light quark content of the nucleon. 

While triangle diagrams can be treated using effective heavy quark nucleon condensate 
instead of performing an explicit one-loop calculation, such a simple treatment is not 
justified in general for box diagrams. Such an approximation would be valid only when 
m q /(Mg — M x ) <C 1 as shown explicitly in [57] for the MSSM. In that case the tree-level 
approach works well for c and b quarks but would fail for t-quarks unless the associated 
squark is much heavier. Nevertheless this is the approximation we use by default in 
micrOMEGAs, the main reason being that in many models the contribution of the Higgs 
exchange diagram is much larger than the one from the box diagrams. The user can 
always ignore this simple treatment and implement a more complete calculation of the 
box diagrams. For example in the case of MSSM-like we have implemented the one-loop 
computation of the neutralino nucleon scattering of Ref. [57], see Section 0] and Appendix 
A. 

In a generic new physics model, new heavy coloured particles can also contribute to 
the WIMP gluon amplitude, for instance squarks in the MSSM. For heavy quarks, the 
computation of the triangle diagrams involving squarks, or any other scalar colour triplet, 
also reduces to a calculation of WIMP-squark scattering with an estimation of the squark 
content in the nucleon. The latter can be obtained by calculating the contribution of 
squarks to the QCD /3-function just as was done for heavy quarks, Eq. HHJ Note however 
that the contribution of scalars to the trace anomaly has an additional factor of 2 due to 
the different dimension of scalar and fermion fields. After substituting A/3 and 7 we get 
at order a s , 

(N\2Ml^ Q \N) = -±- (l + ^ (N\a s G, v G^\N) (48) 

Thus the contribution of scalars is expected to be small because of small scalar content 
in the nucleon. This relation is also known to order a 2 s [§5] . On the other hand other new 
particles such as a heavy Majorana fermion or a real scalar which belong to adjoint color 
representation have very large nucleon densities 

(N\m Q $ Q il> Q \N) = -±-.( N \ aaGlu ,cr\N) (49) 

{N\2MI^q\N) = -±-(N\a s G^GnN) (50) 

In summary, in micrOMEGAs we check the list of coloured particles in the model and 
according to their spin and colour define the nucleon content for each particle using Eq. H7J- 
1501 We then compute the contributions from all WIMP-coloured particles processes. The 
coefficients of the operators for such interactions are calculated automatically in the same 
manner as the coefficients for WIMP-quarks interactions as described in Section 3.1. The 
case of of a color octet vector particle is not treated. 
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4 The special case of supersymmetry 



In the case of the MSSM, additional corrections to both the Higgs exchange diagram and 
the box diagrams with quarks or squark exchange have been computed explicitly. These 
corrections can easily be extended to the case of the CPVMSSM with complex phases 
and of the NMSSM with an extra singlet field. 

The main contribution to the scalar neutralino nucleon cross section involves the Higgs 
exchange diagram with the Higgs coupling to light quark pairs. While in the computation 
of the relic density one could safely neglect the Higgs coupling to light quarks, for direct 
detection this approximation is no longer valid. The MSSM model file must therefore 
be modified accordingly. Furthermore SUSY-QCD corrections to the Higgs exchange 
diagrams can be large and should be taken into account. In particular the gluino-squark 
loops can give large corrections to all Hiqq couplings (if; = h, H, A), they are also related 
to the corrections to the quark masses [78]. These SUSY-QCD corrections will induce 
couplings of the quarks to the 'wrong' Higgs doublet and are enhanced by tan (5 for down 
type quarks. The SUSY-QCD corrections to the Hi — > bb vertices are already taken into 
account in micrOMEGAs for the computation of the relic density of DM and are also 
important for the computation of b — > 57. We generalize this to include the SUSY-QCD 
corrections for all down-type quarks [95]. As in micrOMEGAs_1.3, we define the effective 
Lagrangian 



I m q 1 

C eff = VAna QEDi + A — — 



—Hqq- 



cosa ( _ Am g tana 



cos (3 \ tan (3 

Am, \ If Am n 

+iAqq tan p I 1 — - + hqq 



tan/3 2 / cos/5 V tan a tan /3 



(51) 



where q — d,s, b. We use the same conventions as in [HE] for the supersymmetric model 
parameters and our code is compatible with the SUSY Les Houches Accord (SLHA) [99] . 
For light quarks, Am^ includes only the squark gluino loop while additional electroweak 
corrections which are important for the third generation are included in Amj, [HE] • These 
corrections are much smaller and can be neglected for up type quarks. As we mentioned 
before there is no explicit dependence on the quark mass for the scalar neutralino nucleon 
cross section, nevertheless the effect of SUSY-QCD corrections remains. 

Other vertices that contribute to neutralino quark interactions and also have a depen- 
dence on the light quark mass are the vertices. These vertices involve a coupling of 
the Higgsino component which is proportional to the quark mass as well as a coupling 
of the gaugino component which depends on the squark mixing matrix. Since the off- 
diagonal element of the 2X2 mixing matrix for a given squark flavour is also proportional 
to the quark mass, mixing cannot be neglected even for the first two generations. The 
MSSM model therefore needs to be modified to add these extra terms in vertices for light 
quarks. In general, input parameters for micrOMEGAs are taken from the SLHA [HH] and 
assume no mixing for the first two generations, we define the mixing angle as 

tan20„ = -2 aT9 q « 9 52 
q M? -Ml K ' 

where A u = A u — /i/ tan/3 for up-type quarks and Ad = Ad — /itan/3 for down- type 
quarks, m q is the quark mass including SUSY-QCD corrections [78] • For third generation 
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sfermions, the mixing angle is taken directly from the spectrum calculator as given in the 
SLHA [990. 

The importance of QCD and Arrib corrections in a few sample supersymmetric models 
is displayed In Table El A more detailed description of the parameters of the models 
considered can be found in Tables El El The QCD corrections in Eqs. l4"Tf4*8l are relevant 
only for heavy quarks and will always lead to an increase of the spin independent cross 
section, typically around 5%. The SUSY-QCD corrections on the other hand depend 
on the sign of fi and are enhanced at large values of tan/3. In the decoupling limit, 
these corrections affect mainly the couplings of heavy Higgses as well as the couplings of 
squarks while the light Higgs, which often gives the dominant contribution to the scalar 
cross section remains unchanged. For fi > 0, Am, > and the cross section decreases, 
so that the net effect of the QCD and SUSY-QCD corrections remains small. For fi < 0, 
Am, < and both the heavy Higgs and squark contribution increases. In the sample 
model KP of Table [2] the SUSY QCD corrections lead to a mild decrease of <r„ , this is 
because there is a destructive interference between the dominant light Higgs exchange 
and the heavy Higgs and squark diagrams. 

Table 2: Effect of higher-order corrections on a^ 1 in sample SUSY models. 





BP 


KP 


IP 


NUH 


MSSM1 


tan (3 


10 


40 


35 


30 


10 


fi 


+ 




+ 


+ 


+ 


tth 2 


0.101 


0.101 


0.113 


0.088 


0.100 


a* 1 x 10 V 
tree-level 


8.25 


8.63 


26.1 


9.12 


19.1 


QCD 


8.50 


9.13 


26.8 


9.36 


19.8 


Arrib 


7.53 


8.52 


20.0 


7.75 


18.1 


QCD+Am b 


7.77 


9.02 


20.5 


7.97 


18.7 


QCD+Am 6 +box 


7.78 


9.02 


20.5 


7.97 


18.7 



Finally in the MSSM and its extensions we also compute more precisely one-loop 
corrections. In that case we replace the tree-level treatment of heavy quarks by a one-loop 
computation of the process \g — > X9 including triangle and box diagrams involving quarks 
and squarks. The corrections to the dominant Higgs exchange triangle diagram discussed 
above are always included so the bulk of the corrections are taken into account with 
either option. The one-loop correction from the box diagram is obtained by modifying the 
denominators entering the tree-level diagrams, details can be found in Appendix A. This 
method reproduces the complete results of Drees and Nojiri [57]. We have implemented 
these analytic results and provide a special routine that can be used for comparing the 
two approaches. In the last two lines of Table [21 we compare the results using heavy 
quark currents with the ones including the more complete treatment of loop corrections 
(box diagrams). In mSUGRA models where the Higgs exchange diagrams dominate, 

5 In the case of Isajet, we set Ad = A s = Aj, and A u = A c = A t . Note that this assumption is not 
strictly correct but the mixing in the down-squark sector is in any case dominated by the \i tan (3 term 
while in the up-squark sector it is mostly relevant for the top squark. 
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discrepancies are below the per-cent level. Larger discrepancies can be found in the 
general MSSM, where, based on explicit calculations, we expect loop corrections to be of 
order m q /(Mq — M x ). This can be large for heavy quarks specially when squark masses 
are of the same order as the neutralino masses. 



5 Scattering rates on nuclei: form factors and veloc- 
ity distribution 

To get the rate for direct detection of WIMPs as a function of the recoil energy of the 
nucleus we must take into account both the finite velocity of WIMPs and the nucleus 
form factor. 



5.1 Spin independent interactions 

First consider the simplest case of the SI interaction. We note that ignoring form factor 
effects, we expect isotropic scattering in the center on mass frame. This means that in 
the laboratory frame for a WIMP of velocity v one gets a constant distribution over the 
recoil energy in the interval < E < E max (v). In the non-relativistic approximation, 

E max {v) = 2 f (53) 

For an incoming WIMP with a fixed velocity the recoil energy distribution is thus of the 
form 

dvf si e(E max (v)-E)Fl(q) 

dE ~ a ° E max (v) (54j 

where F A (q) is the nucleus form factor which depends on the transfer momentum q = 

DM particles have a certain velocity distribution, f(v). After integration over incoming 
velocities, the distribution of the number of events over the recoil energy reads 

dJ w = ^w/ 1 ^ {X ' Z + K(A ~ z)f m (55) 

where po is the DM density near the Earth, M^ ef the mass of the detector and t the 
exposure time and 

1(E) = ^-dv (56) 

Jv min {E) v 

EM A ^ 1/2 



v 



v mi n(E) = \-2jf) (57) 



For SI interactions, the form factor is a Fourier transform of the nucleus distribution 
function, 

F A (q) = [ e-^p A (x)d 3 x (58) 
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where p A {x) is normalized such that F A (0) = 1. In our package, we use the Fermi 
distribution function 

/ \ Cnorm /rn\ 

PA{T) = l + exp((r-R A )/a) (59) 

where c norm is fixed by the normalization condition. The resulting form factor is often 
referred to as the Woods-Saxon form factor [33] . The parameters R A and a can be 
extracted from muon scattering data. Compilations of these parameters for several nuclei 
are available [B21E3]- A two-parameter fit to these tables |H] gives 



R A = 1.2343 - 0.6fm (60) 

for a surface thickness, a = 0.52 fm. This is the default value we use for all nuclei. We 
have tested this interpolation. Using the data on the first three radial moments of charge 
distribution [62] we calculate the form factor at a recoil energy of 15 keV, a value slightly 
above threshold for DM detectors. We then calculate the half-density radius R A which 
reproduces this value. Fig. [5] shows that for various nuclei, the value of R A extracted this 
way is well approximated by Eq. 



5.2 Spin dependent interactions 

The SD case features the same kinematics as the case just discussed. The simplest 
approximation would be to assume that the coefficients have the same q-dependence 
as the SI form factor ( 1581) . For a more precise evaluation of the q dependence, one should 
note that when q ^ the vectors J p and J n are not collinear anymore. As a result, 
three form factors need to be introduced. They correspond to the three coefficient of the 
quadratic function of ^ and £ p £ n in the squared amplitude, Eq. [121 Equivalently, one 
can construct the isoscalar and isovector combinations 

ao = C P + £n (61) 
ai = €p- in (62) 

so that the SD recoil energy distribution for a fixed WIMP velocity reads 
da s A D _ 16/4 , a , n 2 , c , , , 2 ME max (v)-E) 



dE 2J A 



-(S 00 (q)a + 5 i(g)a ai + SnO?)^) . (63) 

1 Jl/m.n.m V 



The coefficients Soo(q), Sn(q), and Soo(q) are the nuclear structure functions which take 
into account both the magnitude of the spin in the nucleon and the spatial distribution 
of the spin. They are normalized such that 



S 00 (0) = C(J A )(S p + S n ) 2 

S n (0) = C(J A )(S p -S n ) 2 (64) 
5-01(0) = 2C(J A )(S P + S n )(S p - S n ) 

where C(J A ) = Q±±Mt*±V 

AttJ a 

With this normalization, one recovers the cross section at rest, Eq. [121, starting from 
Eq. 
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After taking into account the velocity distribution, the distribution for the number of 
events over the nuclei recoil energy for spin dependent interactions reads 

^ltE~ = ^Ja + IM ( g °°( g ) a o + 5 oi(g)a ai + S n (q)a 2 1 )I(E) (65) 

The form factors are calculated from detailed nuclear models including the momentum 
dependence. The list of SD form factors implemented in micrOMEGAs is given in the 
Appendix. For nuclei for which the form factor has not been computed precisely, we can 
describe the SD form factors with a Gauss distribution, 

Sij(q) = S t] (0)exp(-q 2 R 2 /4) , (66) 

where «SV, (0) is defined in Eq. [65] and S p and S n can be found in [HI [6]]. To find out 
the effective nucleus radius to be used in the Gauss distribution, we perform a fit to the 
known form factors. First we observe that with a recoil energy E — 15 keV the form 
factors have a similar q dependence (within 10%), 

Sii(g) _ goo(g) _ ggi(g) (67] 
Sn(0) ~ 5 00 (0) ~ Soi(0)' 1 1 

Thus, we can use the same effective radius for all form factors. To extract the A depen- 
dence of R we assume that Z exchange dominates, which means that Sn is the dominant 
form factor. We thus obtain 



R A = 1.7 A 1/3 - 0.28 - 0.78 (A 1/3 - 3.8 + yj (A 1 / 3 - 3.8) 2 + 0.2 ) fm (68) 

A comparison between this approximate formula and the value extracted from the form 
factors tabulated in micrOMEGAs shows that the approximation works quite well, see 
Fig. [2 For light nuclei (A < 100), R A = 1.5A 1/3 fm is also a very good approximate 
formula [67]. Some dependence on the DM model remains because of the different q 
behaviour of Sij. This dependence has been shown to be small [66] . 



Table 3: Number of events per kg/ day x 10 5 for 5 — 50 keV recoil energy for SD interactions 
on various nuclei for model BP. Comparison between two sets of form factors and the 
Gauss approximation in eq. 
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23 Na 
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29 Si 


39 K 


73 Ge 


93 Nb 
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127j 


129 Xe 


131 Xe 
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2.03 
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2.32 
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4.60 


0.79 
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1.55 


Gauss 


16.4 


2.05 


3.36 


2.26 


2.70 


4.67 


2.50 


4.68 


0.75 


6.25 


1.46 


Sxx Tab. 9 




2.13 




2.32 




7.92 




5.54 


1.27 


4.37 


1.37 


Gauss 




2.02 




2.19 




7.81 




6.06 


1.23 


4.60 


1.30 



We compute the number of events for SD interactions for the nuclei listed in Table [3] 
in the region 5 < E < 50 keV. For this we choose the sample model BP in Table [5] and 
compare the results using the precise form factorsof Table 8 and Table 9 [65] with the 
ones obtained by using a Gauss distribution with the approximated formula Eq. [HE] as well 
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A V3 



Figure 2: Ra for SI (full/blue) and SD (dotted/red) Gauss distribution form factors for 
various nuclei, the curves correspond to the approximations [60] and [68j 

as the coefficients S p and S n in Eq. [651 We find that the approximate results with the 
Gauss distribution agree well with the exact form factors, the agreement is usually better 
than between different sets of form factors. Even for small nuclei where the interpolation 
for the radius Ra does not work as well there is a very good agreement between the 
different sets of form factors, see Table El This is because for light nuclei the momentum 
is very small hence the form factor does not need to be known precisely. Note that the 
Gauss approximation would not work as well at large recoil energies where the form factor 
decreases too rapidly. 

5.3 Velocity distribution of dark matter 

The nuclear recoil energy measured in direct detection experiments depends on the WIMP 
velocity distribution in the rest frame of the detector (I55|I65I) . This in turn depends on 
the WIMP velocity distribution in the rest frame of the galaxy and the Earth velocity 
with respect to this frame. The latter is determined from various observations. The 
measurements of the velocity of Sun and others objects close to the Sun give a value both 
for the velocity of rotation of the LSR(local standard of rest) [HE] 

v = 220 ± 20km/ s (69) 

and for the peculiar velocity of Sun in this system PI EO] E 

v pec = (10.0, 5.2, 7.2)km/s (70) 

6 Here we use the galactic co-ordinates (X,Y,Z) where X is toward the Galactic center, Y in the direction 
of rotation and Z toward the north Galactic pole. 
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The Earth velocity with respect to the galactic frame is thus the sum of v = (0, v , 0), v pec 
and of the Earth velocity in the solar system. Assuming that the Earth's orbit is circular 
and that the axis of the ecliptic lies in the Y-Z plane, the Earth velocity in Galactic 
co-ordinates is 

v e = v e ( — sin(27r£), sin7cos(27rt), cos7cos(27rt)) (71) 

where v e = 29.79km/ s, 7 = 30.5° and t is in years. More precise expressions, taking into 
account an elliptic orbit and the orientation of the axis of the ecliptic can be found in 

The velocity distribution of DM particles on the Earth is obtained from the DM 
velocity distribution in the rest frame of the Galaxy, Fqrf, 

f(v) = J 5{v - \V\)F GRF (V -v - v pec - v e )d 3 V (72) 

Because the mass of Galaxy is finite there is some v max such that Fq RF = for \ V\ > v max , 
astronomical observations [72j give the 90% confidence interval 

498km/ s < v max < 608km/ s 

with a median likelihood of v max = 544km/ s. 

There are several models of DM velocity distribution [73], they are correlated with the 
DM density distribution. The simplest and most widely used model to describe the DM 
density is the isothermal sphere model [H]. In such a model the DM velocity distribu- 
tion corresponds to a Maxwellian distribution. In our package we have implemented a 
truncated Maxwellian distribution, 

F G rf(V) ~ exp(-\V\ 2 /AV 2 )Q(v max - \V\) (73) 

which leads to 



vi) 2 \ ( miniv + v u v max ) 2 



ex P aTTo — — ex P 



AV 2 J "V AV 2 
where AV = Vq in the isothermal model, c norm is fixed by the normalization condition 

f(v)dv = 1 



(74) 



and 

vi = \vq + v pec + v e \ w v + (v pec ) y + v e sin 7 cos(27rt) 

Note that the Earth motion around the Sun leads to a 7% modulation effect of V\ and 
in turn to a modulation of the signal in direct detection experiments. This modulation 
was investigated in [7TI [73] for a large set of models of DM spatial distributions. The 
DM velocity distribution close to the Sun could be quite different from the Maxwell 
distribution. For example condensation of cold DM in clumps and streams will lead 
to a delta-function distribution. This function has also been implemented in the code. 
The impact of the different velocity distributions on the limit on the detection rate of 
neutralinos was also discussed in [71] . In order to allow for deviations from the isothermal 
model, in the implementation of the Maxwellian distribution in micrOMEGAs we treat AV, 
Po, Vi as free parameters. Furthermore alternative models for DM distribution can always 
be implemented by the user. 

The DM density near the SUN is estimated to be in the range po — 0.1—0.7 GeV/cm 3 |75j . 
As default, we take the commonly used value po = 0.3 GeV/cm 3 . 
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6 Description of routines 



All routines are available both in C and Fortran and in both cases they have the same 
names and the same arguments. For simplicity, here we describe only C routines. The 
types of parameters are listed in 

sources/mi cromegas.h for C routines and in sources/mi cromegas_f.h for Fortran routines. 

• setProtonFF (scalar , ps_vector, sigma) 

• setNeutronFF (scalar , ps_vector, sigma) 

These two routines allow to redefine the default values of the coefficients describing 
the quark contents of nucleons, see Section [331 Each parameter is an array of type double 
of dimension 3 whose elements correspond to the form factors for d, u, and s quarks 
respectively. Default values of form factors are displayed in Table HI To keep the default 
value for one set of coefficients, use NULL (NOFF in Fortran). 

Table 4: Default values of quark form factors in proton and neutron. 



current 


proton 


neutron 


d 


u 


s 


d 


u 


s 


scalar 


0.033 


0.023 


0.26 


0.042 


0.018 


0.26 




-0.427 


0.842 


-0.085 


0.842 


-0.427 


-0.085 




-0.23 


0.84 


-0.046 


0.84 


-0.23 


-0.046 



• getScalarFF (m u /mrf, m s /md, ct^n, 0"o,FF_proton,FF_neutron ) 

This routine computes the scalar coefficients for quark contents in the nucleon from the 
mass ratios m^/m^, m s /mrf as well as from a n N an d c"o following Eq. I2TH2TI1 It can be used 
to find the input values of the quark coefficients in setProtonFF and setNeutronFF. a n ^ 
and do are specified in MeV. The coefficients are stored in FF_proton and in FF_neutron, 
both arrays of dimension 3 and of type double. 

• FeScLoop(sgn, mq,msq,mne) 

This function computes the loop K-factor for MSSM-like models, Eq. IA-51 with a spin 
1/2 WIMP and a scalar "squark". It allows to replace the denominators for s,w-channel 
diagrams as specified in eq. IA-41 Here mq, msq,mne refer to quark, squark, and WIMP 
masses. 

• nucleonAmplitudes(qBOX,pAsi ,pAsd,nAsi ,nAsd) 

This routine calculates amplitudes for WIMP-nucleon elastic scattering at zero momen- 
tum. pAsi (nAsi) are spin independent amplitudes for protons (neutrons), whereas pAsd(nAsd) 
are the corresponding spin dependent amplitudes. Each of these four parameters is an 
array of type double and dimension 2. The first (zeroth) element of these arrays gives 
the x- nuc l eon amplitudes whereas the second element gives x-nucleon amplitudes. Am- 
plitudes are normalized such that the total cross section for either \ or x cross sections 

a tot = ^(\A SI \ 2 + 3\A SD \ 2 ) (75) 

7T 

The qBOX parameter specifies a function like FeScLoopO that is used to improve 
the tree-level calculation in order to reproduce the results of the box-diagram calculation. 
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The function FeScLoop included in our package can be used for spin 1/2 WIMPS and 
scalar "squarks". To obtain the tree-level result, substitute qBOX = NULL (qBOX = 
NoLoop) in C(Fortran). nucleonAmplitudes returns a value different from zero only 
when there is an internal problem in calculation. 

• MSSMDDt est (loop, &pS , &pV , &nA , &nV) 

This routine computes the proton (neutron) scalar, pS(nS) and vector, pV,(nV) one- loop 
or tree-level amplitudes in the MSSM using directly the formulae of Drees and Nojiri [57] . 
If loop = then calculations are done at tree level assuming heavy quark condensates 
in the nucleon, otherwise the one-loop results for neutralino interactions with gluons are 
used. The QCD and SUSY-QCD corrections are included. Amplitudes are normalized 
according to ( 1751) . This routine exists only in the C- version. 

• fDvMaxwell(v) 

returns f(v)/v where f(v) is defined in Eq. ( 1741 . The argument v is expressed in km/s. 
This function is used as argument of the nucleusRecoil function described below. 

• Setf Maxwell (DV,vl , vmax) 

sets parameters for f DvMaxwell. The arguments correspond to the parameters A V, v±, 
and v max in Eq. (1741) in km/s units. Unless this routine is called the program uses the 
default values (220,225.2,700). 

• SetfDelts(v) 

sets velocity of DM for ^-function distribution. 

• fDvDelts(v) 

indicates to nucleusRecoil that the velocity distribution is a delta function which is non- 
zero for the parameter specified in Setf Delts (v) . 

• SetFermi(C,B,a) sets parameters for A-dependence of the Fermi half radius, Ra = 
CA 1 ' 3 + B and for the surface thickness, a. Default values are given in Eq. [601 

• nucleusRecoil (rho , f Dv , A , Z , J , S00 , SOI , Sll , qBOX , dNdE) 

This is the main routine of the direct detection module. The input parameters rho, f Dv, 
specify the DM velocity distribution while A, Z, J, S00, SOI, Sll specify properties of 
detector material. The return value gives the number of events per day and per kilogram 
of detector material. The distribution over recoil energy is stored in the array dNdE. 
This array has to be of type double with 200 elements. The value in the i th element 
corresponds to 

dN. 

, n E=i*keV 

dE 

in units of (1/keV/kg/day). For a complex WIMP, nucleusRecoil averages over WIMP 
and WIMP. 

The input parameters are: 
rho - density of DM near the Earth in GeV/cm 3 ; 

f Dv - this parameter is a function which specifies the DM velocity distribution, f Dv= 
f{v)/v where f(v) is defined in Eq. (172|) . the velocity v is given in km/s and fDv(v) in 
(s/km) 2 . The function f DvMaxwell described above is the default function that can be 
used here. 

A - Atomic number of nucleus; 

Z - Number of protons in the nucleus; 

J - nucleus spin. 

qBOX - a parameter needed by nucleonAmplitudes, see the description above. 

S00(p), S01(p), Sll(p) are nucleus form factors for spin-dependent interactions. They 
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are functions of the momentum transfer in fm~ l (argument of double type). These form 
factors are assumed to be normalized as in Eq. [65l 

Nucleus form factors for a wide set of nuclei are implemented in micrOMEGAs. The 
available form factors are listed in Tables El El in Appendix B and are extracted from the 
review article of Bednyakov and Simkovic [HH]. These functions are named 

S{xx}{Nucleus Name}{Atomic Number}[A] (76) 

where xx is either 00,11 or 01. The last character is optional and is used to distinguish 
different implementations of form factor for the same isotope. 

For convenience our package has predefined constants for nucleus charge 

Z_{Name} 

and nucleus spin 

J_{ Name } { atomic -number } 
for all isotopes listed in Tab l8|9l 

For example, a call to nucleusRecoil for 73 Ge detector should be 

N=nucleusRecoil (0.3, f DvMaxwell , 73 , Z_Ge , J_Ge73 , S00Ge73 , S01Ge73 , SllGe73 , FeScLoop , dNdE) ; 
• int PlotSS(Sxx,A, title, Emax) 

allows a visual check of the spin dependent form factors. This routine opens a new 
window where Sxx is displayed as function of the recoil energy. Here Sxx stands for any 
of functions in Tab. [H [9J A is the atomic number, title defines the text that will be 
displayed on the screen and Emax is the upper limit of the recoil energy for the plot. 

The PlotSS routine is based on CalcHEP [76] plot facilities. By a click of the mouse 
on one point in the plot area one gets information about the value of X/Y coordinates 
for that point. To close the window press the 'Esc' key. Pressing some other key opens 
the menu which displays min and max values of the function and allows to change limits 
of the Y axis and switch between linear/logarithmic scale of Y axis. For additional help 
press the Fl key. 

•nucleusRecoilO (0.3, f Dv , A , Z , J , Sp , Sn , qBOX , dNdE) 

is similar to the function nucleusRecoil except that it uses Eq. IM1 to define 500(0), 
511(0), 501(0) and the q-dependence of these form factors is the one associated with the 
Fermi distribution with radius, Ra given in Eq. [UH This function can for example be 
used for light nuclei, such as 3 He, or for any other nuclei where the more precise form 
factor has not been included in Tab. El For all nuclei listed in Tab. El as well as for 1 H 
(S p = 0.5, 5 n = 0) and 3 ife [100j . form factors have been predefined with names 

Sp{Nucleus Name}{Atomic Number} and Sn{Nucleus Name}{Atomic Number} 

(77) 

These constants correspond to the form factors listed in Tab JH] extracted from the review 
article [HI] • They all satisfy EqJMl Table El shows how well the approximation works. 

Two auxiliary routines are provided to work with the energy spectrum computed with 
nucleusRecoil and nucleusRecoilO. 
•displayRecoilPlot (dNdE, title, El, E2) 

plots the generated energy distribution dNdE. Here title is a character string specifying 
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the title of the plot and E1,E2 are minimal and maximal values for the displayed energy. 
E1,E2 have to be an integer value in keV units, E2 has to be smaller than 200. This 
routine has the same service facilities as PlotSS described above. Note that the plots can 
be saved in Latex format as well as in GNUPLOT or PAW formats. 
•cutRecoilResult(dNdE, El, E2) 

calculates the number of events in an energy interval defined by values E1,E2 in keV 
units. 

7 Sample results 



Table 5: SI and SD cross sections in sample SUGRA models. 
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X? 


248.4 


97.9 


161.7 


207.4 


141.3 


223.5 


476.3 


109.4 


218.6 


xt 


466.8 


183.0 


302.6 


395.7 


264.4 


284.1 


893.3 


156.7 


404.3 


k 


257.1 


107.6 


170.5 


213.8 


152.3 


2129. 


816.7 


1607. 


254.4 


h 


952.9 


362.5 


649.9 


774.6 


581.8 


1730. 


1836. 


999.3 


848.2 


h 


111.5 


111.0 


112.5 


114.2 


112.4 


118.6 


119.2 


113.7 


115.3 


A 


888.2 


420.5 


576.4 


769.0 


426.2 


1547.6 


927.4 


1614. 


497.8 


Qh 2 


0.126 


0.101 


0.115 


0.107 


0.113 


0.101 


0.128 


0.119 


0.088 


a bI x 10 9 pb 
set A 


0.98 


7.50 


2.48 


0.019 


19.6 


9.03 


0.509 


46.9 


7.65 


set B 


1.82 


15.9 


4.98 


0.055 


44.1 


14.4 


1.06 


85.7 


16.7 


set C 


0.77 


5.49 


1.87 


0.011 


13.9 


7.67 


0.375 


37.2 


5.50 


a b v u x 10 6 pb 
set A' 


0.301 


5.00 


1.97 


0.360 


3.78 


140. 


0.096 


510. 


2.22 


set B' 


0.230 


3.69 


1.56 


0.261 


3.05 


129. 


0.082 


470. 


1.89 


a b n u x 10 6 pb 
set A' 


0.251 


4.34 


1.55 


0.322 


2.89 


87.0 


0.069 


317. 


1.58 


set B' 


0.334 


5.89 


2.01 


0.443 


3.69 


99.1 


0.085 


362. 


1.93 


N xl0 3 /kg/day 
73 Ge 


0.19 


3.00 


0.70 


0.0071 


6.08 


2.41 


0.056 


20.3 


1.68 


131 Xe 


0.31 


5.35 


1.18 


0.0085 


10.5 


3.36 


0.089 


31.6 


2.79 



The predictions for the spin independent cross sections on protons and spin dependent 
cross sections on neutrons and protons are listed in Table [5] in the case of the MSSM with 
input parameters fixed at the GUT scale. The first seven models (AP to MP) are inspired 
by the mSUGRA benchmarks of Ref . [101] . The input parameters have been modified such 
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that the prediction for Qh 2 falls near the central WMAP value when the top quark mass is 
fixed to the value measured at Tevatron, m t = 171.4 GeV |102j . The parameters of Model 
BP were adjusted to the ones of the SPA1A benchmark |103j . The spectrum calculator 
used is SuSPECT |104] . The last two models are non-universal SUGRA models. In model 
NUG, the gaugino masses are non- universal at the GUT scale with M 3 = 0.7M 2 = 0.7 Mi 
leading to a lighter coloured sector than in the universal case. In model NUH only the 
Higgs masses are non-universal with Mh u = 2mo and Mu d = — 0.6mo- For each model 
the predictions for the SI cross section are given for three sets of values for the coefficients 
describing the quark density contents in the nucleon. Set A correspond to the default 
values, Eq. [35], while the input parameters for set B (a w N = 70, er = 35 MeV) and set 
C (o- wN = 55, a = 40 MeV) are varied within the expected range, Eq. [2S Note that 
these do not correspond to the most extreme choices, yet the predictions for the cross 
sections can vary by a factor of 2 or 3. For the SD cross sections, set A' and set B' 
correspond respectively to the old (Eq. 1391) and new (Eq. 00]) estimates of the quark 
density coefficients. In Table |5] the predictions for the number of events per day and 
per kg of detector material for 73 Ge and l31 Xe are also presented. Here we assume the 
default values for the quark density coefficients as well as for the velocity distribution. 
The distribution of the number of events as function of the energy is also computed. The 
results for models BP and KP for a detector made of 131 Xe are displayed in Fig. [3] 




20 40 60 80 100 120 140 160 180 200 



E(keV) 

Figure 3: dN/dE for a 131 Xe detector in two mSUGRA models specified in Table El BP 
(full) and KP (dash). 

We have also compared our results with other public codes. For this we have used 
the same spectrum calculator (here we have taken Isajet_7.75) and have removed from 
micrOMEGAs the QCD and SUSY QCD corrections which are neglected in other codes. 
For a given set of quark coefficients, our results for a SI are in very good agreement with 
Isajet_7.75 for the dominant contribution due to Higgs exchange, differences appear in the 
squark exchange contribution. For SUGRA models this can lead to 25% differences be- 
tween the two codes. For a SD , our prediction is usually below the one of Isajet, differences 
between the two codes can reach a factor 4. This can be traced back to a sign discrepancy 
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between the Z and the squark exchange diagram 0. As concerns DarkSUSY_4.1, we have 
discrepancies that can reach 25% for a SI , this is due to the Higgs exchange diagram. In 
DarkSUSY the running mass is used in the hqq vertex while the pole mass is used when 
estimating the quark scalar coefficient, whereas in micrOMEGAs we use the same mass 
everywhere as we have discussed in Section 3.3. Our predictions for a SD are in very good 
agreement with those of DarkSUSY for the Z exchange contribution. We however have a 
factor 2 difference in the squark exchange diagram which can lead to 50% discrepancies 
in a SD for our SUGRA test models. Increasing the squark exchange contribution by a 
factor 2 in DarkSUSY, we recover very good agreement with micrOMEGAs 0. 



Table 6: Comparison with Isajet_7.75 and DarkSUSY_4.1 in sample mSUGRA models 
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In the sample models of Table [7] the dependence on the axial vector coefficients for the 
spin dependent cross-section is far less important (within 30% between sets A' and B') 
than for the scalar cross section. Furthermore much smaller variations are observed if one 
just sticks to the uncertainty associated with the latest experimental results, Eq. HUJ This 
is specific to SUGRA models where one finds that the Z exchange diagram completely 
dominates over the squark exchange diagram. This is valid in all models where Z exchange 

7 The direct detection module is being improved in Isajet, we expect much better agreement between 
the two codes in the next public version, Isajet_7.76 |105j . 

8 The DarkSUSY code has been updated and is now in good agreement with micrOMEGAs [106j 
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Table 7: SI and SD cross sections on protons in non-SUGRA models, all masses in GeV. 
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dominates. In the more general MSSM with weak scale input parameters, it is possible 
to find models where the squark exchange diagrams contribute significantly to the spin 
dependent cross section and can furthermore interfere with the Z exchange diagram. This 
can lead to a strong reduction of either the proton or neutron cross section as shown for a 
sample MSSM model in Table [7] where the most relevant parameters are specified in the 
Table, in addition we set m & = M A = -A t = l TeV, m~ lh = 500 GeV and M 3 = 800 GeV. 
Finally sample results for the case of a gauge boson in the Little Higgs model (LH) and 
of a Dirac right-handed neutrino (RHN1 and RHN2) are listed in Table [71 Note that in 
this table the values for SI cross sections on neutrons only are given, in most models these 
numbers are similar to the ones of protons except in the RH neutrino model where the 
cross section on protons is very much suppressed. In this model the small dependence 
on the quark coefficients is due to the Higgs exchange diagram. The results for the RH 
neutrino model agree with the ones of 



3 We use the version of the Little Higgs model implemented into CalcHEP and available at 



http: / /hep. pa. msu.edu more details can be found in |107j . 
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8 Installation 



The package can be obtained from the web page wwwlapp . in2p3 . f r/lapth/micromegas. 
All instructions for installing the package can be found in |108j . Contrary to pre- 
vious versions, the user must download only one file micromegas_2.2.tgz which con- 
tains the full implementation of different models such as the MSSM, CPV-MSSM |l(J9j . 
NMSSM [IIQ], [III] as well as the little Higgs model (LHM) [W] and the right-handed 
neutrino model(RHNM) [38]. Facilities to incorporate other new models remain as in 
micrOMEGAs_2.0, specific instructions on how to install new models are given in [108J. 
Note that the identifiers SO and V5 are now reserved for the auxiliary fields used to gen- 
erate automatically the operators in Table [H they cannot be used to designate particles 
in the model. Note also that since micrOMEGAs searches for quarks in the list of par- 
ticles by identifying the PDG codes, care has to be taken to implement these correctly 
in the model file. The direct detection module is not compatible with earlier versions 
of the supersymmetric models, since light quark masses and the corresponding couplings 
of light quarks to Higgses need to be introduced. To run the code we provide for each 
model one sample file, main, c (main.F in Fortran). This sample program can be used to 
compute the relic density of DM, the cross sections for direct and indirect detection, the 
cross sections at colliders and decay widths. Various options can be set in that program 
depending on the need of the user, the various switches available are commented in the 
main file. The sample programs cycle2.c and cycle5.c found in the MSSM directory will 
reproduce the numerical results in Tables I2ll5l 

9 Summary 

The new module for computation of the cross-section for WIMP scattering on nucleus in 
micrOMEGAs described here applies to any type of CDM candidate, whether Majorana 
or Dirac fermion, scalar or vector boson. After the new model is implemented within 
micrOMEGAs, the mass, the spin and the interactions of the WIMP candidate are computed 
automatically. Because the nucleon has an important light quarks component, WIMP 
interactions with light quarks in the nuclei often dominate over those of heavy quarks. 
For scalar interactions the amplitude for WIMP quark scattering is proportional to the 
quark mass, thus the mass of light quarks have to be taken into account and incorporated 
into the model file. In the specific case of the MSSM this means expanding the model 
file since for the relic density calculation all fermions of the first two generations were 
taken to be massless. Furthermore within this model, higher-order corrections to the Hqq 
vertices for down-type quarks such as SUSY-QCD corrections are taken into account. For 
all models, we use an effective vertex for WIMP interactions with heavy quark in the 
nuclei, this takes into account dominant QCD corrections to Higgs exchange diagrams. A 
more precise treatment is available for MSSM-like models. 

We have also explicitly shown on a few examples the impact of varying the coefficients 
for the quark density content of the nucleon. Uncertainties are very large for scalar 
interactions and are in general more under control for spin dependent interactions that 
arise through Z exchange. For other contributions to spin dependent interactions the 
uncertainties can be quite large although the value predicted is often several orders of 
magnitude below the present limit. 
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Note added 

An online tool for computation of direct detection rates with micrOMEGAs 2.2 has been 
set up by Rachid Lemrani, see 

http : //pisrvO .pit .physik.uni-tuebingen.de/darkmatter/micromegas_g/. 

Appendix A - Box diagrams 

As shown explicitly by Dress and Nojiri, the tree-level calculation overestimates s-channel 
amplitudes when M q — M x < m q . In that case one should instead compute completely 
the box diagrams or use the procedure described here. In the MSSM, the Lagrangian for 
neutralino-quark-squark interaction has the generic form 

C Y = q{a + bj 5 )xq + h.c. (A-l) 

and leads to the WIMP-quark scattering amplitude at tree-level 

I . b 2 a 2 > 

~ 4 Wf - [M x + m q y M? - (M x - m q f> ^ ' > 

The explicit computation of the box diagram [57J leads to an amplitude of the form 

A box = l - ((b 2 - a 2 )F D (m q , M q , M x ) + (a 2 + b 2 )F s (m q , M q , M x )) (A-3) 

where Fp and F$ are loop functions. A simple modification of the propagators in IA-21 
will therefore reproduce the amplitude A box , 

K{±l,m q ,M^M x ) (A-4) 



M 2 -(M x ±m q ) 2 
where 

K(s, m q , Mq, M x ) = F D (m q , M q , M x ) + sF s (m q , M q , M x ) 

= ^(m/i-^-sM^-^-^)) (A-5) 
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and 



" x x 2 -2x + 2/3 .. e . 

h = / dx Tyj — — ( A "6) 



D 2 



r 1 , x(x 2 - 2x + 2/3) 
h = / cfe^ '-L (A-7) 



Is = (A-8) 
J 4 = / V 3(1 D ~ X)2 (A-9) 



1 Jm x(1 -x){2-x) 

Of 1 



h = I ^ (A-10) 



D = x 2 Ml+x{M 2 ~-m 2 q - M 2 )+m 2 q (A-ll) 

The modification of the denominators of tree level diagrams in micrOMEGAs is achieved 
by a small modification of the CalcHEP code. Obviously this trick will work in any MSSM- 
like model where the WIMP is a fermion that interacts with quarks via scalar quarks, in 
particular in the NMSSM and the CPV-MSSM. The same trick can be generalized to other 
models, for this however one has to compute the appropriate loop factors which depend 
on the WIMP and 'squark' spin. Note that the box diagrams for WIMP scattering on 
gluons includes also diagrams where gluons couple to squarks. When the WIMP scattering 
amplitude is computed with the loop K-factor option, micrOMEGAs omits the tree- level 
processes involving squarks (or in general colored triplet that are not quarks) assuming 
that their contributions are included in the K-factor. 



Appendix B - Spin dependent nucleus form factors 



Table 8: Nucleus SD form factors realized in micrOMEGAs 



Identifier 


Isotope 


Ref 


data in [65] 


comments 


SxxF19 


19 p 




Eq. 7 






SxxSi29 


29 St 


m 


Eq. 14 






SxxNa23 


23 Na 


[66 


Eq. 9 






SxxTel25 


125 Te 


[66] 


Eq. 18, Tab. 


IV 


Bonn-A potential 


SxxI127 


127j 


[66] 


Eq. 20 




Bonn-A potential 


SxxXel29 


129 Xe 


[BB] 


Eq. 21 Tab. 


IX 


Bonn-A potential 


SxxXel31 


131 Xe 


[BB] 


Eq. 21 Tab. 


IX 


Bonn-A potential 


SxxA127 


27 Al 


H13] 


Eq. 11 






SxxK39 






Eq. 15 






SxxGe73 


73 Ge 




Eq. 17 






SxxNb92 


93 Nb 


|115| 


Tab. II 




Scanned plot 
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Table 9: Alternative nucleus SD form factors. 



Identifier 


Isotope 


Ref 


data in [65J 


comments 


SxxSi29A 

SxxNa23A 

SxxTel25A 

SxxI127A 

SxxXel29A 

SxxXel31A 

SxxGe73A 

SxxXel31B 


29 St 
23 Na 

125 Te 
127j 

129 Xe 
131 Xe 
73 Ge 
131 Xe 


Uifil 

1112J 

P] 

P] 

[66] 

[66] 

|116| 

|117| 


Eq. 12 
Eq. 10 

Eq. 18, Tab.IV 
Eq. 19, Tab.VI 
Eq. 21 Tab. IX 
Eq. 21 Tab. IX 
Eq. 16 
TABLE VII 


Nijmegen II potential 
Nijmegen II potential 
Nijmegen II potential 
Nijmegen II potential 

5 1 31 (0) should be 0.04^ 



Appendix C - Additional routines 

To facilitate model independent comparisons with data, we provide additional routines to 
compute the nucleus recoil energy using as input the WIMP mass and the cross sections 
for SI and SD scattering on nucleons. 

• nucleusRecoilAux(rho , f Dv, A, Z, J, S00 ,S01 , Sll ,Mwimp, cslp, csln, csDp, csDn,dNdE) 
This function is similar to nucleusRecoil and returns the number of events per day 
and per kg of detector material. The additional input parameters include the WIMP 
mass, cslp (csln) the SI cross sections for WIMP scattering on protons (neutrons) and 
csDp(csDn) the SD cross sections for protons (neutrons). A negative value for one of 
these cross sections corresponds to a destructive interference between proton and neutron 
amplitudes. 

• nucleusRecoilOAux (rho , f Dv , A , Z , J , Sp , Sn , Mwimp , cs Ip , csln , csDp , csDn , dNdE) 
This function is similar nucleusRecoilAux except that it uses Eq. [64] to define 500(0), 
511(0), 501(0). For details see the description of nucleusRecoilO. 

• setRecoilEnergyGrid(step,nStep) 

This function redefines the grid for the recoil energy. After calling this function the 
nucleusRecoil function returns a recoil energy distribution where = i x step keV, 
i — 0, nStep — 1. 
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